This directory has programs that compute 2 and 3 point correlation functions.
They are hosted at 

    http://code.google.com/p/mjarvis/.

If you have any problems or discover a bug, please post an issue about it 
there.  If you have other questions about using the code, you can also feel 
free to email me at:

    michael@jarvis.net.

Note: if you use this code for research that makes it into a paper, the 
appropriate reference is: 

    Jarvis, Bernstein, & Jain, 2004, MNRAS, 352, 338.  

The code is in transition from an old version where most parameters were hard
coded, so you were required to recompile whenever you changed anything to a 
new version where parameters are specified in a parameter file.  The two-point 
code has been converted, so this Read.me file is about that.  However, I have
not yet converted the three-point code to this system, so if you want to run 
that, please see Read.me.old which describes the old programs.


********
Overview
********

The program that computes the two-point correlation function is called corr2.
It includes several varieties of correlations:
NN = the normal two point correlation function of things like 2dF that
     correlate the galaxy counts at each position.
NE = correlation of counts with shear.  This is what is often called
     galaxy-galaxy lensing.
EE = two-point shear correlation function.

There are plans to add scalar values as well, which we plan to call K for 
kappa, since weak lensing is my focus, but it could also be used for 
temperature values or any other scalar quantity.  If you would like to use 
this code for that application, please post an issue to help push this
closer to the top of my TODO list.


************
Installation
************

Currently you need to manually edit the makefile to set your compiler, 
directories, etc.  It should be fairly straightforward.  

The primary items you will want to set are:

CC = The C++ compiler to use.  Default = g++
OMP = The appropriate OpenMP flag for that compiler.  Default = -fopenmp
      To disable OpenMP support, you may comment out this line entirely.
FITSDIR = The directory where cfitsio is installed.  Default = /usr/local
BINDIR = The directory where you want the corr2 executable installed.
         Default = .

There are some other items below that which you probably won't need to edit,
but might if you need to include different library names or something strange
like that.

Then, to build the code, just type:

make

As usual, you can do make clean to remove all the .o files to save disk space.


*************
Running corr2
*************

The executable corr2 takes one required command-line argument, which is the 
name of a configuration file:

corr2 config_file

A sample configuration file is provided, called default.params, and see below 
for a fairly complete list of 
includes much of this 

You can also specify parameters on the command line after the name of 
the configuration file. e.g.:

corr2 config_file file_name=file1.dat e2_file_name=file1.out
corr2 config_file file_name=file2.dat e2_file_name=file2.out
...

This can be useful when running the program from a script for lots of input 
files.


*******************************
Parameters about the input file
*******************************

fie_name = (str or list) The file(s) with the galaxy data.

    You an also specify two files here, in which case the program calculates a 
    cross-correlation between the two sets of values.  e.g.
    file_name = file1.dat file2.dat

    If you are specifying this on the command line, you'll need to put 
    quotes around the names, or it won't be parsed correctly:
    filename="file1.dat file2.dat"


do_auto_corr = (bool, default=false) Whether to do auto-correlations within
               a list of files.
do_cross_corr = (bool, default=true)) Whether to do cross-correlations within 
                a list of files.

    If there are more than two names in the file_name paramter, then the code 
    will normally calculate all the pair-wise cross-correlations of each pair 
    i,j in the list.  The default is to only do cross-correlations, but not 
    auto-correlations.  You may change that behavior by changing do_auto_corr 
    or do_cross_corr.  


file_name2 = (str or list) The file(s) to use for the second field for a 
             cross-correlation.

    If you want to cross-correlate one set of files with another, then the 
    above file_name parameter isn't sufficient.  Instead, you would list 
    the first set of files in file_name, and the second set in file_name2.
    Of course, each list may only contain one file each, so this is another
    way to specify two files to be cross-correlated.


file_type = (ASCII or FITS) The file type of the input files.

    The default file type is normally ASCII.  However, if the file name 
    includes ".fit" in it, then a fits binary table is assumed.
    You can override this behavior using file_type.


first_row = (int, default=1)
last_row = (int, default=-1)

    You can optionally not use all the rows in the input file.
    You may specify first_row, last_row, or both to limit the rows being used.
    The rows are numbered starting with 1.  If last_row is not positive, it 
    means to use all the rows (starting with first_row).

x_col = (int/str) Which column to use for x
y_col = (int/str) Which column to use for y
ra_col = (int/str) Which column to use for ra
dec_col = (int/str) Which column to use for dec

    For the positions of the objects, you can specify either x,y values, which 
    imply a flat-sky approximation has already been performed (or ignored),
    or ra,dec values, which are of course positions on the curved sky.

    For ASCII files, the columns are specified by number, starting with 1 being
    the first column (not 0!).  
    For FITS files, the columns are specified by name, not number.

x_units = (str, default=arcsec) The units of x values.
y_units = (str, default=arcsec) The units of y values.
ra_units = (str) The units of ra values.
dec_units = (str) The units of dec values.

    All distances on the sky include a "units" parameter to specify what in 
    units the values are specified.  Options for units are radians, hours, 
    degrees, arcmin, arcsec.  If omitted, arcsec is assumed if you are using 
    x,y.  But for ra, dec the units field is required.


g1_col = (int/str) Which column to use for g1
g2_col = (int/str) Which column to use for g2
e1_col = (int/str) Which column to use for e1
e2_col = (int/str) Which column to use for e2

    If you are doing one of the shear correlation functiond (i.e. NE or EE),
    then you need to specify the shear estimates of the corresponding galaxies.
    The g1,g2 values are taken to be reduced shear values.  If your input 
    catalog lists distortions rather than shears, the program can convert
    for you.  Just indicate this by declaring columns are e1_col and e2_col.

w_col = (int/str) Which column to use for the weight (if any)

    The weight column is optional. If omitted, all weights are taken to be 1.

flip_g1 = (bool, default=false) Whether to flip the sign of g1
flip_g2 = (bool, default=false) Whether to flip the sign of g2

    Sometimes there are issues with the sign conventions of gamma.  If you 
    need to flip the sign of g1 or g2, you may do that with flip_g1 or flip_g2 
    (or both).

Note: If you are cross-correlating two files with different formats, you may 
      set any of the above items from file_type to flip_g2 as a two element 
      list.  In this case, the first item refers to the file(s) in file_name,
      and the second item refers to the file(s) in files_name2. However, you 
      may not mix (x,y) with (ra,dec) or (g1,g2) with (e1,e2).  The former 
      because it wouldn't make sense.  The latter just because it's not 
      enabled in the code.


project = (bool, default=false) Whether to do a tangent plane project for
    handling the curved sky values.

    The native RA, Dec code is almost as fast as the flat-sky code, so it is
    generally preferable to let the code handle RA and Dec using the correct
    spherical geometry formulae.  But if you want, you can instead project the 
    RA, Dec values onto a tangent plane.  This is probable most useful for 
    investigating bugs in the curved sky code, so I suspect most users will
    not want to use this feature.

project_ra = (float) The ra of the tangent point for projection.
project_dec = (float) The dec of the tangent point for projection.

    The default tangent point for the projection is the average position.
    However, this may be inappropriate for some reason, so you may specify the 
    projection point with project_ra, project_dec.
    (These use the same units specified for ra_units, dec_units.)


*****************************************************************
Parameters about the binned correlation function to be calculated
*****************************************************************

min_sep = (float) The minimum separation to include in the output.
max_sep = (float) The maximum separation to include in the output.
nbins = (int) The number of output bins to use.
bin_size = (float) The size of the output bins in log(sep).

    The bins for the histogram may be defined by setting any 3 of the above 4 
    parameters.  The fourth one is automatically calculated from the values
    of the other three.

    There is one exception.  If you set min_sep, max_sep, and bin_size, 
    then it won't generally be the case that the corresponding number of 
    bins is an integer.  So the code will increase max_sep slightly to make
    sure the total range is an integer number of bins.

sep_units = (float, default=arcsec) The units to use for min_sep and max_sep.

    sep_units is also the units of R in the output file.

bin_slop = (float, default=1)

    The code normally determines when to stop traversing the tree when all of the 
    distance pairs for the two nodes have a spread in distance that is less than the 
    bin size.  i.e. the error in the tree traversal is less than the uncertainty 
    induced by just binning the results into a histogram.  This factor can be changed
    by the parameter bin_slop.  It is probably best to keep it at 1, but if you want to
    make the code more conservative, you can decrease it, in which case the error 
    from using the tree nodes will be less than the error in the histogram binning.
    (In practice, if you are going to do this, you are probably better off just 
    decreasing the bin_size instead and leaving bin_slop=1.)

    Note, if you set bin_slop=0, then the code will effectively do a brute-force
    calculation, since it will branch all the way to each leaf of the tree.

smooth_scale = (float)

    In addition to the raw output, the code will also optionally output a smoothed 
    version of the correlation functions, which is better for plotting.
    The smoothing scale is specified as smooth_scale.
    If omitted or smooth_scale = 0, then no smoothing will be done.


***********************************
Parameters about the output file(s)
***********************************

The kind of correlation function that the code will calculate is based on 
which output file(s) you specify.  It will do the calculation(s) relevant for 
each output file you set.  For each output file, the first line of the output 
says what the columns are.  See the descriptions below for more information
about the output columns.

CAVEAT: The error estimates for all quantities only include the propagation 
        of the shot noise and shape noise through the calculation.  It 
        does not include sample variance, which is almost always important.
        So the error values should always be treated as an underestimate
        of the true error bars.

n2_file_name = (str) The output filename for point-point correlation function.

    This is the normal density two-point correlation function.  The output 
    columns are:

    R_nominal = The center of the bin

    <R> = The mean separation of the points that went into the bin.  
          Technically, since we bin in log(R), this is really exp( <log(R)> ).

    omega = The NN correlation function = (DD-DR-RD+RR) / RR

    sig_omega = The 1-sigma error bar for omega

    DD,DR,RD<RR = The 4 direct numbers from which omega was calculated.
                  Note: For an auto-correlation, DR and RD are identical, but 
                  for cross-correlations, they may be different.


ne_file_name = (str) The output filename for point-shear correlation function.

    This is the point-shear correlation function, often called galaxy-galaxy
    lensing.  The output columns are:

    R_nominal = The center of the bin

    <R> = The mean separation of the points that went into the bin.  
          Technically, since we bin in log(R), this is really exp( <log(R)> ).

    <gamT> = The mean tangential shear with respect to the point in question.

    <gamX> = The shear component 45 degrees from the tangential direction.

    sig = The 1-sigma error bar for <gamT> and <gamX>.

    weight = The total weight of the pairs in each bin.

    npairs = The total number of pairs in each bin.

    R_sm (if smooth_scale is set) = <R> for the smoothed values

    gamT_sm (if smooth_scale is set) = <gamT> smoothed over the appropriate scale.
    
    sig_sm (if smooth_scale is set) = The 1-sigma error bar for gamT_sm.


e2_file_name = (str) The output filename for shear-shear correlation function.

    This is the shear-shear correlation function, used for cosmic shear.
    The output columns are:

    R_nominal = The center of the bin

    <R> = The mean separation of the points that went into the bin.  
          Technically, since we bin in log(R), this is really exp( <log(R)> ).

    xi+ = <g1 g1 + g2 g2> where g1 and g2 are measured with respect to the
          line joining the two galaxies.

    xi- = <g1 g1 - g2 g2> where g1 and g2 are measured with respect to the
          line joining the two galaxies.

    xi+_im = <g2 g1 - g1 g2>.  In the formulation of xi+ using complex 
             numbers, this is the imaginary component. 
             It should normally be consistent with zero, especially for an
             auto-correlation, because if every pair were counted twice to 
             get each galaxy in both positions, then this would come out 
             exactly zero.

    xi-_im = <g2 g1 + g1 g2>.  In the formulation of xi- using complex 
             numbers, this is the imaginary component.
             It should be consistent with zero for parity invariant shear 
             fields.

    sig_xi = The 1-sigma error bar for xi+ and xi-.

    weight = The total weight of the pairs in each bin.

    npairs = The total number of pairs in each bin.

    R_sm (if smooth_scale is set) = <R> for the smoothed values

    xi+_sm, xi-_sm (if smooth_scale is set) = xi+, xi- smoothed over the 
                   appropriate scale.
    
    sig_sm (if smooth_scale is set) = The 1-sigma error bar for xi+_sm, xi-_sm.


m2_file_name = (str) The output filename for the aperture mass statistics.

    This file outputs the aperture mass variance and related quantities, 
    derived from the shear-shear correlation function.

    Note: for the definition of the aperture mass, we currently use the 
    Crittenden choices for U, Q:
      
        U = 1/2Pi (1-r^2/2) exp(-r^2/2)
        Q = 1/4Pi r^2 exp(-r^2/2)

    To switch to the Scneider choices:

        U = 9/Pi (1-r^2) (1/3-r^2)
        Q = 6/Pi r^2 (1-r^2)
    
    you can #define SCH_U in the file CalcT.cpp and recompile.  I may
    eventually switch to having this be specifiable via a configuration 
    parameter.  But I haven't done so yet.

    The output columns are:

    R = The radius of the aperture.  (Spaced the same way as  R_nominal is 
        in the correlation function output files.

    <Map^2> = The E-mode aperture mass variance for each radius R.

    <Mx^2> = The B-mode aperture mass variance.

    <MMx>(a), <MMx>(b) = Two semi-independent estimate for the E-B cross term.
                         (Both should be consistent with zero for parity
                         invariance shear fields.)

    sig_map = The 1-sigma error bar for these values.

    <Gam^2> = The variance of the top-hat weighted mean shear in apertures of 
              the given radius R.

    sig_gam = The 1-sigma error bar for <Gam^2>.


nm_file_name = (str) The output filename for <N Map> and related values.

    This file outputs the correlation of the aperture mass with the 
    aperture-smoothed density field, derived from the point-shear correlation 
    function.

    R = The radius of the aperture.  (Spaced the same way as  R_nominal is 
        in the correlation function output files.

    <NMap> = The E-mode aperture mass correlated with the density smoothed
             with the same aperture profile as the aperture mass statistic
             uses.

    <NMx> = The corresponding B-mode statistic.

    sig_nmap = The 1-sigma error bar for these values.


norm_file_name = (str) The output filename for <N Map>^2/<N^2><Map^2>
                       and related values.

    This file outputs the <N Map> values normalized by <N^2><Map^2>.  This 
    provides an estimate of the correlation coefficient, r.
    
    R = The radius of the aperture.  (Spaced the same way as  R_nominal is 
        in the correlation function output files.

    <NMap> = The E-mode aperture mass correlated with the density smoothed
             with the same aperture profile as the aperture mass statistic
             uses.

    <NMx> = The corresponding B-mode statistic.

    sig_nm = The 1-sigma error bar for these values.

    <N^2> = The variance of the aperture-weighted galaxy density.

    sig_nm = The 1-sigma error bar for <N^2>.

    <Map^2> = The aperture mass variance.

    sig_mm = The 1-sigma error bar for <Map^2>.

    nmnorm = <NMap>^2 / (<N^2> <Map^2> )

    sig_nmnorm = The 1-sigma error bar for this value.

    nnnorm = <NN> / <Map^2> 

    sig_nnnorm = The 1-sigma error bar for this value.


************************
Miscellaneous parameters
************************

verbose = (int, default=0) How verbose the code should be during processing.

    0 = no output
    1 = normal output
    2 = extra output


num_threads = (int, default=auto) How many (OpenMP) threads should be used.

    The default is to let OpenMP determine an appropriate number of threads 
    automatically.  Usually this matches the number of cores your system has.


**************
Reporting bugs
**************

If you find a bug running the code, please report it at:

    http://code.google.com/p/mjarvis/issues/list

Click "New Issue", which will open up a form for you to fill in with the
details of the problem you are having.

If you would like to request a new feature, then after clicking "New Issue",
you can change the Template to "Feature Request", which will ask more
appropriate questions about the feature you would like me to add.

